Correlation of spin and velocity in granular gases 
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In a granular gas of rough particles the spin of a grain is correlated with its linear velocity. We 
develop an analytical theory to account for these correlations and compare its predictions to numer- 
ical simulations, using Direct Simulation Monte Carlo as well as Molecular Dynamics. The system 
is shown to relax from an arbitrary initial state to a quasi-stationary state, which is character- 
ized by time-independent, finite correlations of spin and linear velocity. The latter are analysed 
systematically for a wide range of system parameters, including the coefficients of tangential and 
normal restitution as well as the moment of inertia of the particles. For most parameter values 
the axis of rotation and the direction of linear momentum are perpendicular like in a sliced tennis 
ball, while parallel orientation, like in a rifled bullet, occurs only for a small range of parameters. 
The limit of smooth spheres is singular: any arbitrarily small roughness unavoidably causes signif- 
icant translation-rotation correlations, whereas for perfectly smooth spheres the rotational degrees 
of freedom are completely decoupled from the dynamic evolution of the gas. 

PACS numbers: 45.70.-n, 45.70.Qj, 47.20.-k 



I. INTRODUCTION 

Materials which are composed of macroscopic objects, 
i.e. granular media, attract increasing scientific inter- 
est due to their importance in nature and technology, 
e-g- [11, 01 ■ The latter may be exemphfied by trans- 
port and storage of sand, cereals, granular chemicals, 
etc. the former-by avalanches, land slides, dust devils, 
etc. Spectacular celestial objects, like planetary rings 
or interstellar dust clouds, can serve as another exam- 
ple of natural granular systems The granular mat- 
ter exists there in a gaseous state and exhibits many 
properties of a common molecular gas, e.g. P, 0, H, H]. 
The main (and very important) difference of a granu- 
lar gas from a molecular gas is the dissipative nature of 
particle interactions, which are describe by macroscopic 
mechanics of solids rather than by a microscopic inter- 
action potential. The consequences of the dissipative 
interactions are quite substantial: A spatially homoge- 
neous state is unstable 0, H, velocities are not dis- 
tributed according to a Maxwell-Boltzmann distribution 
[13, M, [H, [H, M, [H, [H, p, [li and the diffusion or 
self-diffusion is anomalous |l9l. |20|. [2ll . [2^ . [23j . These 
properties of a granular gas have been observed for the 
case of smooth particles, when grain collisions do not 
affect their rotational motion. This is, certainly, an over- 
simplified model, since real grains have a rough surface 
and exchange rotational and translational energy in col- 
lisions. 

Real granular particles experience frictional forces 
when colliding. Hence, a more adequate model takes 
into account the rotational motion of particles and the ex- 
changeof rotational and translational en ergy in collisions 

[3, iM 113, m, m, 113, m, m, no, im, m, mm, hi . Dis- 
sipative frictional gases exhibit additional unusual fea- 



tures which arc not present in molecular gases. For 
instance, equipartition between rotational and transla- 
tional motion does not hold [13] and the hydrodynamic 
description requires an additional field and a dynamic 
equation to account for its evolution [1^, [3^. Moreover, 
the rotational and translational motion of particles are 
correlated as mentioned in a very implicit way in Ap- 
pendix E of [26j and worked out in [33] • In the present 
study we analyze the latter effect in detail. 

In Sec [n] we introduce a model of frictional particles 
and the observables of interest. Subsequently in Sec. 
mil an approximate analytical theory is developed and 
in Sec. IIVI we briefly explain the simulation techniques. 
The main results are presented in Sec. |Vl where we com- 
pare predictions of the analytical theory with data from 
simulations. The emphasis lies on the correlations in the 
quasi-steady state, but we also briefly discuss the relax- 
ation to the steady state. The technical details of the 
calculations are given in the Appendix. 



II. MODEL AND OBSERVABLES 

We consider a granular gas consisting of inelastic 
hard spheres of radius a, mass m, and moment of inertia 
/ = qma^. Here the dimensionlcss variable q is deter- 
mined by the mass distribution within the disc. The 
state of the system is fully described by the particles' po- 
sitions {vi}, velocities {f ;}, ^iid angular velocities {uJi} 
for i = 1, . . . , A. The particles move freely in between 
instantaneous collisions, whereupon their linear and an- 
gular velocities change according to the collision rule: the 
relative velocity at the point of contact of colliding par- 
ticles is 



Vi - V2 + aii X {ivi + 1^2) 



(1) 



2 



with n = ni2 = (ri — / |ri — r2|. The post- 
colUsional (primed) velocity is related to the pre- 
coUisional one by 



g ■n= e„ g ■ n 
g' X fi = Et g X n . 



(2) 



The cocfEcient of normal restitution is denoted by £„ 
with < £„ < 1. The value £„ = implies no rela- 
tive motion in the normal direction after the collision, 
whereas for e„ = 1 no dissipation of the normal com- 
ponent of the relative motion occurs. The coefficient 
of tangential restitution has two elastic limits, namely 
£4 = 1 corresponding to smooth spheres and £* = — 1 cor- 
responding to perfectly rough (reflecting) collisions with- 
out loss of energy for the tangential motion. For all other 
values energy is lost in the tangential component. In gen- 
eral, both coefficients of restitution, £„ and St, depend on 
the impact velocity [sSj, i3^, i4l| . 

Together with the conservation of linear and angular 
momentum the collision rule, Eq. ([2]), determines the 
post-collisional velocities in terms of the pre-coUisional 
ones: 



v[ — Vi — S, u)[ = u)i -\ {n X S) 

qa 

V2 =V2 + S, = u;2 + — (n X ^) , 

qa 



(3) 



where mS denotes the exchange of linear momentum with 



^ = Vtg + {Vn - ■nt){n-g)fi, 



Vn 



g 1 - 



'^*=2 1 



(4) 



(5) 



In the present study we address only non-driven sys- 
tems. Moreover, we focus on the homogeneous cooling 
state (HCS) of a gas, which is characterized by two time- 
dependent granular temperatures, one for the transla- 
tional and one for the rotational motion, 



T 
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m \ - 

?— 1 



and R 
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(6) 



One generally observes that after a transient period 
the system reaches a quasi-stationary state where r = 
R{t)/T{t) = const., that is, both temperatures decay 
with the same rate. In general, r 7^ 1 so that equiparti- 
tion is violated. The value of r depends on the collision 
parameters as well as on the moment of inertia [24l . [4^ . 

In this paper we focus on the correlation between the 
axis of rotation of a granular particle and the direction of 
its linear velocity, which may be quantified by the angle 
9i between the linear and rotational velocity, 



cos 9i 



All information on the angle is contained in the distribu- 
tion 



1 ^ 

/(cos0) = — 5{cos6 — cos6'i). 



(8) 



In a molecular gas all values of cos 9 occur with equal 
probability due to equipartition. In contrast for a gran- 
ular gas we know that equipartition is violated and we 
expect to observe deviations from the equi-distribution. 

Because of symmetry, the average of cos 9i over all par- 
ticles vanishes. Thus, a measure of correlations is the 
second moment, 



(9) 



If the angular and linear velocities are not correlated in 
their direction, (cos^ 6*) = 1/3. Hence, any deviation 
of (cos^ 6*) from 1/3 indicates correlations. Moreover, 
if (cos^ 6*) < 1/3 the angular and linear velocities are 
preferably perpendicular, like in a sliced tennis ball, while 
for {cos^9) > 1/3 they are preferably aligned like in a 
rifled bullet. 



III. ANALYTICAL THEORY 

The evolution of any observable 

Fit)^Fi{u{t),v,it),uj,it)}) (10) 

may be obtained by means of the pseudo-Liouville oper- 
ator via 



dtF{t) = iL+F{t) for t > 0. 



(11) 



For hard spheres the pseudo-Liouville operator decom- 
poses into two parts, = £jq + Ju'_^, where iLg = 
£jq'^ + £,0°' describes the free streaming of translational 
and rotational motion of particles. Here = Vi ■ Vi 
and a similar expression for £;q°* . The latter is not needed 
here, because we never specify the orientation of our par- 
ticles , which are perfect spheres. The interaction part 
of the pseudo-Liouville operator reads, = J2i<ri '^jj^ 
where the binary collision operator 7ij reads [ill . [24 | 

iJij = ~hij ■ VijQ (-n,j • Vij) 5 {r^j - 2a) (bij - ij . 

(12) 

The operator 6^ replaces unprimed by primed values ac- 
cording to the collision rule, Eq. For example, 

bi2Vi = v[ , 6i2f2 = V2 , bi2Vk = ffc , k^l,2 

(13) 

with v[ and V2 given by Eq. Q and with similar relations 
for the rotational velocities. 

The ensemble average of a dynamic variable is defined 

by 

{F), = / drp{0)F{t) = / dTp{t)F{0) (14) 



3 



with dT = l\^(d^r,(Pv,d^uJi). Here F{t) = 
exp {—iL+t)F{0) and p{t) = exp (— liL^^t)p(O) denotes 
the A'^-particle distribution, whose evolution is governed 
by the adjoint ZL^ of the evolution operator Differ- 
entiating Eq. (fTi|) one obtains 



d_ 

di 



{F)t = / dTp{0)j^F{t) = I dTp{Q)ii:+F{t) 
= j dTp{Q)eyi]){iL+t)iL+F{Q) 
= j dTp{t)iSi+F{Q) = {i^i+F)^ . (15) 



It is impossible to compute the time-dependent A^- 
particle distribution exactly, so that we have to resort 
to approximations. A standard procedure in the analyti- 
cal treatment of granular gases is to assume homogeneity 
and molecular chaos, e.g. (see also Under these 

assumptions the A'^-particle velocity distribution function 
takes the form 



=ffJv(ri,...,rjv)]^pi(fj,Wj,i), 



(16) 



where the iV-particle correlation function of a hard 
sphere system, 5Ar(ri, . . . , tat), is not affected by the par- 
ticle roughness. For the HCS it may be approximated by 
the corresponding function of an equilibrium hard-sphere 
system (e.g. Q). For an isotropic system pi{v,(jj) de- 
pends in general on v = \v\,uj = |w| and the angle 9 
{cos 9 = V ■ u)/ {\v\ \uj\). Here we are particularly in- 
terested in the dependence on cos 9 and expand pi in 
Legendre polynomials P„ (cos 6) 



pi(.,u;,t)occxp(-^j cxp(^-^j 



^&„(i)t;"a;"P„(cos0), 



(17) 



n=0 



where the 6„ {t) are time dependent expansion coefficients 
and the distribution function has to be normalized ac- 
cording to J dT pi = 1. We use a simple Gaussian even 
though the distributions are non-Gaussian for strong dis- 
sipation and high densities. Deviations have been han- 
dled by an expansion in Soninc polynomials [iij . Here 
we concentrate on the dependence on cos 9 and leave a 
more general ansatz with both, angular correlations and 
non-Gaussian distributions, to future work. To keep the 
calculations tractable, we limit the calculation to the low- 
est non-trivial order 



pi{v,u},t) (xexp 



exp 



2i?(t)y (18) 



2T{t) ) 

where hit) = h^it) and P2(cos6') ^ 3/2(cos2 9 - i/a). The 
terms for odd n vanish by symmetry. 



The lowest order coefficient h(€) is simply related to the 
quantity of interest (cos^S)^. Using Po(cos0) = 1 and 
expressing (cos^ 9) ^ in terms of Legendre polynomials we 
can write 



(cos^ 9)^^^JJ [Pa{cos9) + 2P2(cos0)] 



(19) 



where for brevity we introduce the shorthand notation 



f— ) 

\2ttTJ 



3/2 



d V exp 



2T 



(20) 



and similarly for . The angular integration in the Eq. 
(|19p may be performed using the orthogonality relation 
for Legendre polynomials, yielding 



( 



cos 



6T{t)R{t) 



(21) 



Hence, the correlations of interest manifest themselves 
through the coefficient 6(t)-the larger the coefficient, 
the more pronounced are deviations from the value 
^cos^ 6*) = 1/3 of the uncorrelated case. 

To summarize our analytical approach so far: The time 
dependent A^-particle distribution has been parametrised 
by three time-dependent functions T{t), R{t) and b{t), 
which have to be calculated self-consistently. This is 
achieved by applying the general equation (|15p for the 
evolution of an observable to T(t), R{t) and h{t) and 
using our ansatz for p{t), see Eqs. (|16ll8p . Even with 
all these simplifying assumptions, the analytical calcula- 
tions are rather cumbersome and all the details of the 
calculation have been relegated to the Appendix. 

The results are three first order differential equations 
for T(t),R{t) and b{t). These simplify, if we measure 
times in units of the Enskog collision frequency loe = 
16(7rT/m)^/^7ia^g2(2a). In other words we rescale time 
according to dr = ujEdt and obtain: 



B 



b{T) T{t)R{t) 



qm?a? 



dR 
d7 



BT{t) - C 



b{T) T{t)R{t) 



qm a'^ 



R{t) 
R{r) 



(22) 



where 



A = r]n{l - rjn) + rjt{l 
Vt 

q 



Vt) 



B 



9 

vt 



c 



(23) 



and 



9 2 

qm a 



^(i) + 5(i) 



T{r) 
40 dR 



40 dT 
T(j}'dT 



ir) 



nT)R{T) 



^(0) + ^(0) 



R{r) 
T{r) 



^(0) 



T{t) 
R{t) 



(24) 
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FIG. 1: (Color online) The coefficients in Eq. (|25|) as a func- 
tion of £t for £„ = 0.9 and q = 2/5. Note that -S'^' is slightly 
negative for Et > 0.7. Except for the coefficient A^"' all coef- 
ficients vanish in the limit et ~* 1- 



The constants are given by: 
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- 1 77„(7;„ - 1) (25a) 



(MD) calculations to check the predictions of the 
analytical theory. DSMC determines the stationary dis- 
tribution of the scaled velocities by numerically solving 
the kinetic Boltzmann equation which is based on the as- 
sumption of molecular chaos. Consequently, for its appli- 
cation it is assumed that the gas is uniform, thus, spatial 
correlations of the particles are neglected. If this pre- 
condition is given, DSMC yields very precise statistical 
results because of the large number of particles which can 
be simulated (here we use TV = 2 x 10'' particles [ioj). 

Molecular Dynamics calculates the trajectories of the 
particles using the collision rule, Eq. ([2]), therefore, MD 
allows to trace the evolution of the correlation. On the 
other hand, MD is restricted to much smaller systems 
as compared to DSMC. Although MD is free from the 
mentioned assumptions, DSMC is significantly more effi- 
cient for a homogeneous granular gas. Moreover, in the 
limit of low density both methods provide, in principle, 
identical results for the stationary state [i^ . In practice, 
we use MD for N = 8000 particles to study the transient 
process of the system's relaxation to its steady-state and 
up to iV = 10^ for steady state correlations. The volume 



fraction is 



N 47ra-' 



V 3 

gas is always in the HCS 



— 0.0146 or even smaller, such that the 
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(25b) 
(25c) 



+ 22ir,t + rjn)~6iVt+vl) (25d) 



V. RESULTS 

Starting from a random distribution of velocities and 
angular velocities with mean (v) = (uj) = 0, after some 
transient period the system relaxes to a steady state 
where the correlation of the spin and the translational ve- 
locity as well as the ratio of translational and rotational 
temperatures adopt stationary values. We quantify these 
correlations by means of the second moment (cos^ 0), see 
Eq. ^ and analyze this quantity as a function of three 
parameters, £„, St, and q in Sec. IV Al The relaxation to 
the steady state is discussed in Sec. IV Bl and in Sec. IV CI 
we consider correlations beyond the second moment and 
investigate the distribution of cos 6*. 
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3 q 
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(25e) 



Eqs. and constitute a set of self-consistent equa- 
tions for the observables T{t),R{t), and b{t). Fig. [1] il- 
lustrates the dependence of the above coefficients on the 
coefficient of tangential restitution. 



IV. SIMULATIONS 

We performed both Direct Simulation Monte Carlo 
(DSMC) [1^ and event-driven Molecular Dynamics 



A. Steady-state correlations 

To study the steady-state properties it is convenient to 
introduce an auxiliary variable 



xir) 



qinn?a? 



(26) 



Using x{t) and r(T) = R{t)/T{t) we recast the set of 
three equations (|22I24[) for b, R and T into a set of two 
equations for x and r. The result reads 



5 



dr 
dx 



B-C 



x{t) 



r{T) + Ar(r) - B 



x{t) 



r2(r) 



-x{t) 1^(1' + B(^V(t) - 2QA - 20C + 20B[r(r) + r-i(r)] + 202;(t)[C - Br(T)]/2| 
A(0)-S(")r(r)-C("V-i(r). 



(27) 
(28) 



Setting the left hand side of Eqs. ^\ and to zero 
one arrives at a set of coupled nonlinear equations for the 
stationary values Too = r{T oo) and Xoo = x{t — > oo). 
Instead of solving these equations directly, we resort to 
an iteration scheme: At the outset we calculate a first 
approximation of the temperature ratio neglecting 
correlations, that is, for x = 0. Hence we assume that 
for moderate inelasticity and roughness the tempera- 
ture ratio is not noticeably affected by the rotational- 
translational coupling. The result reads 



0.335 



.(0) 



A~C 
2B 



Using this value for the stationary temperature ratio we 
then proceed to calculate an approximate value of Xoo 

A(i)+B(i)rS^-40C + 40B/r^^ 
where we use the fact, that 

=^-C + 2i3/r(^) (31) 

and neglect the terms quadratic in Xoo since they are 
presumably small. In principle, one could further iterate 
to get better approximations, but we find that the results 
are reasonably good already at this stage. For the more 
intuitive variable, cos^ 9, Eq. ([501) implies 



(cos^ e)^ 



1 6 +5WrS^ 
3 5^(1) +40B/r. 



(0) 



(0) 

OO 



40C 
(32) 

Fig. [2] shows the steady-state value of the correlation 
factor (cos^ 6)^ as a function of £t for different values of 
En in comparison with DSMC results. Obviously, theory 
as well as simulations show that both types of correla- 
tions may occur, (cos^ 6*) < 1/3, as for a sliced tennis ball 
or (cos^ 6*) > 1/3 as for a rifled bullet The dependence of 
the correlations on St is nonmonotonic with the strongest 
correlations for et ^ and et — ^ 1. Even though the de- 
pendence on Sn is also not strictly monotonic, the domi- 
nant tendency is an increase of correlations with decreas- 
ing e„, i.e. increasing inelasticity. The agreement be- 
tween theory and computer experiment is excellent for 
small inelasticity. Moreover, even for significant dissi- 
pation the theory is able to reproduce qualitatively the 
simulation results. 
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FIG. 2: (Color online) Steady-state value of (cos^ ^)oc ^ 
function of the coefficient of tangential restitution, St, for dif- 
ferent e„. The predictions of the analytical theory, Eq. (|32p . 
are depicted by lines and points indicate the simulation data 
by DSMC. The line of vanishing correlations, (cos^ 0) — 1/3 
is shown, as well as the isolated point et = 1, which refers to 
the system of perfectly smooth hard spheres. Note the exis- 
tence of non- vanishing correlations even in the limit of smooth 
spheres, ej — > 1 (see Eq. I33|l . 



Decreasing the moment of inertia, q, turns the mag- 
nitude of the correlations more sensitive to changes in 
the coefficients of tangential restitution, as one can see 
from Fig. [31 Interestingly, varying the moment of inertia 
can even alter the type of the correlations: For instance, 
for (7 = 1/5 there exists a region for ej > 0, where the 
rotation axis is preferably directed along the linear veloc- 
ity, (cos^ 9)^ > 1/3, while for q = 2/3 there is no such 
region. 

Fig. m (upper panel) illustrates the analytical result, 
Eq. ([32l) for the whole range of parameters St and £„. 
Note that for the majority of values of the coefficients, 
(cos^ 6*)^ < 1/3, that is, in most cases the axes tend to be 
perpendicular to each other. Only in two small regions 
of the parameter space the axes are preferably parallel. 
The correlations vanish only for combinations of e„ and 
Et indicated by full lines. Dashed lines show curves of 
constant r. Strong correlations appear for large devia- 
tions from equipartition. This is shown more clearly in 
the middle and bottom panels of Fig. [H which demon- 
strate the rather strong influence of the moment of inertia 
/ on the correlation factor (cos^ 9) . 
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FIG. 3: Steady-state value of (cos^ 6^)^ for e„ — 0.9 as a 
function of £t and for different moments of inertia of a grain 
(see also Fig. |8]). With the decreasing moment of inertia 
the correlations become more sensitive to variations of the 
coefficient of tangential restitution 



To check the assumption that strong correlations occur 
for strong deviations from equipartition, we plot in Fig. [5] 
the correlation factor (cos^ Q) as a function of T^a and 
£t. Technically this may be done, using e„ = e„(roo)- 
thc inverse function of r^a = ?'oo(en), given by Eq. (j29p . 
for each fixed St- Note that pronounced correlations are 
present mainly for strong dissipation and large temper- 
ature ratios. Also note the small range of admissible 
temperature ratios for very rough spheres. 

Analyzing Eq. (j32p in the limit of vanishing roughness, 



et-l^ 3 8 7 



(33) 



we see that even the smallest roughness induces finite 
correlations, for any given (fixed) value of the coefficient 
of normal restitution, £„ 7^ 1. For = 1, that is, for 
perfectly smooth spheres, the initial rotational velocity 
of the particles is preserved. Therefore, the initial rota- 
tional energy is preserved as well and r does not reach a 
steady state. On the other hand (cos-^ 6*) relaxes to the 
stationary value 1/3 once the correlations in the initial 
values of the translational velocities are lost due to colli- 
sions. Hence a straightforward expansion around et = 1 
is problematic, or at least should be done with much 
care, as long as there is a finite inelasticity e„ 1. [See 
also the discussion of relaxation times in the following 
paragraph.] 



B. Relaxation to the steady-state 

So far we have discussed the quasi-stationary state, 
which is characterized by constant r and (cos^6'). It is 
also of interest to understand, how this stationary state 
is reached — starting from arbitrary initial conditions. 
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FIG. 4: Stationary value (cos^ 6)^ (color coded) as a function 
of normal (£„) and tangential (st) coefficients of restitution. 
The stationary value of the temperature ratio r is superim- 
posed through the dashed contour lines. The solid lines in- 
dicate vanishing correlations {{cos 9}^ = 1/3). The moment 
of inertia q = 2/5 (upper panel) corresponds to homogeneous 
spheres. The middle and bottom panel show the same data 
for g = 1/5 and q — 2/3, respectively 



Of particular interest is the limit of almost smooth 
spheres 77^ oc et — 1 ^ 1 [see the definition. Eq. ([5])]. 
While the decay of the rotational temperature R and 
the translational temperature T takes place extremely 
slowly, that is, with a rate ~ 774 <C 1 [see Eq. ([^ with 
r = R/T ~ A/B in this limit], the relaxation of the 
temperature ratio, r = R/T as well as of the correlation 
factor X ~ ((cos-^ 6') — 1/3) occurs on the collision time 
scale. Indeed, in this limit one can write using Eqs. (j27p . 
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FIG. 5: Stationary value of (cos^ O) ^ (color coded) as a func- 
tion of the temperature ratio r^o and the coefficient of tan- 
gential restitution et- As previously, the solid lines indicate 
vanishing correlations and the dashed lines follow constant 
values of the coefficient of normal restitution £„. Note the 
logarithmic scale for the r-axis. The ragged border is an ar- 
tifact of the limited numerical resolution 



and the definitions of the coefficients (^5]) (^5]) . 
dr/dr c -ry„(l - 77™) ('Vr£)) (r - r^) , 

with r^^ ~ A/B - I/t^I > 1 from Eq. This implies 

that r relaxes to its stationary value exponentially fast 
with a rate ?7n(l — Vn) = 0(1) (that is, on the collision 
time scale), while both temperatures T and R continue 
to decay with the same small rate. 

To analyse the relaxation of x{t) to its steady state 
value, we use Eq. (pS)) and approximate r(T) by its steady 
state value : 



dx 
d7 



-ao — aix — a2X 



where 

ao = 
ai ^ 



1 

20 
1 

20 

1 



1 



(34) 



(35) 



^(1) 



20 



A-C + B{r^ 



02 = -(C-Broo). 



The above equation with the initial condition a;(0) = is 
solved by 



1 — tanh ( 
with the relaxation time 



1 — tanh 



Trol 



Trcl 



4aoa2 



(36) 



(37) 



and tanh(/) = ai/ \/a\ — 4aoa2. Evaluating the coeffi- 
cients for typical values of et and e„, we find that the 




FIG. 6: Relaxation time Ty-^i (in the collision units) 
of ^cos^ when it approaches the steady-state value 
(cos^ '^ote the narrow range of possible values for Trd 



relaxation of the correlation factor <^cos^ 0)^ to its steady- 
state also occurs within a few collisions per particle. This 
is illustrated in Fig. [6l where we plot the relaxation time 
Troi given by Eq. 

We wish to stress here again, that the relaxation on 
the collisional time scale to the steady state values ap- 
plies only to the temperature ratio and the mean square 
cosine of the angle between linear and angular velocity. 
For nearly smooth particles, £t — > 1, the relaxation of 
the rotational and translational temperatures is, never- 
theless, a very slow process, which proceeds with a small 
rate, tending to zero as et ^ 1. 

To demonstrate the existence of several time regimes 
we discuss in the following an instructive example. We 
initialize the particles with ui = corresponding to r = 0. 
The collision parameters are e„ = £t = 0.8 so that the 
asymptotic value of the ratio of temperatures is Too > 
1. We expect r to monotonically increase as a function 
of time — and this is indeed observed as shown in Fig. 
[71 Now, we can check our hypothesis that correlations 
are small for values of r close to equipartition. If the 
hypothesis is correct, we should observe non-monotonic 
behavior of (cos^^)^. For short times the correlations 
should be large and of tennis ball type, because grazing 
collisions are the most effective for spinless particles to 
gain angular momentum. At intermediate times, when 
r I, the correlations should be very small or vanishing. 
In the asymptotic state with Too > 1, one should again 
observe finite correlations. 

These three time regimes are clearly born out in the 
time dependent correlations, shown in Fig. [Tj (a) In the 
short time regime (0 < t < lO'^) correlations are strong 
and < r < 1. (b) At intermediate times (10^ <t< 10^) 
equipartition holds approximately r sa 1 and correlations 
are small or vanishing, (c) The steady state {t > 10^) 
is characterized by r ^ 1 and finite (cos^S)^ < 1/3. 
The agreement between analytical theory and molecular 
dynamics is good also for the time-dependent quantities. 
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0.28 





FIG. 7: (Color online) Relaxation of (cos^ 6)^ and of the ratio 
of temperatures r{t) = R(t)/T(t) to the steady state. Dots: 
molecular dynamics data for 8000 particles, lines: analytical 
theory. To show that vanishing correlations (cos^ 6)^ coincide 
with equipartition, we have chosen the vertical axes, such 
that the point r = 1 on the right ax;is (blue) and the point 
{cos^e)^ = 1/3 on the left axis (red) have the same vertical 
height as indicated by a horizontal line. 



FIG. 9: The ratio of the temperatures calculated with h{t) 
according to Eq. (|24p to those with h(t) = 0. The coefficients 
of restitution are £„ = St = 0.8, and the initial ratio of ro- 
tational to translational temperatures was set to the steady- 
state value r(0) — r^o- The inset shows r(t)/r(0) — 1 as a 
function of time. Note that the deviation of r{t) from roa is 
always very small. 



0.335 




t 



FIG. 8: Impact of the grains' moments of inertia on the evo- 
lution and steady-state of (cos^^)^. The system parameters 
are Sn = 0.9, et ~ 0.9, r(0) = 0.001 and vanishing initial 
correlations. The values oi q {I — qma?) represent spheres 
with the mass concentrated towards the center (g = 1/5), the 
homogeneous spheres {q = 2/5) and spheres with the mass 
concentrated mainly in the outer shell {q — 2/3). 



Fig. [5] demonstrates that the moment of inertia of 
the partieles does not change the evolution of (cos^ 6*)^ 
quahtatively. For the particular choice of the coefficients 
of restitution the correlations are more pronounced for 
larger q ~ I /ma? and fade with decreasing q. This how- 
ever is not a general rule; depending on the coefficients 
e„ and et, this tendency may reverse. 

The correlations between translational and rotational 
motion also have a noticeable, albeit small impact on the 



basic characteristics of granular gases — the translational 
and rotational temperatures. In Fig. [5] we present the 
time dependence of R{t)/R'^^\t) — the ratio of the rota- 
tional temperature R{t) with correlations to the corre- 
sponding value i?^*'^ (t) without correlations. The respec- 
tive ratio T{t)/T^^'>{t) for the translational temperature 
is also plotted. Here we choose the ease of large r ~ 24 
(e„ = Et = 0.8), which correspond to (cos^S)^ < 1/3, 
that is, for preferably perpendicular rotational and trans- 
lational velocity. Fig. [5] demonstrates that the effect 
of the correlations on the granular temperatures R{t) 
and T(t) is indeed small. The corresponding quantity 
r{t) = R{t)/T{t) is also not sensitive to these correla- 
tions. Moreover r{t) does not deviate noticeably from 
its steady-state value throughout the system's evolution, 
that is, |r(r)/T'oo — 1| ^ 1, as shown in the inset ofFig.[9l 



C. Beyond the second moment 

A complete one-particle picture includes the distribu- 
tion 

1 ^ 

W(cos 0, v,u}) — — ^ (5(cos 9 — cos Oi)S{v ~ Vi)d{Lu — uji). 

m 

Since correlations are developed in collisions, one intu- 
itively expects that particles with larger velocities, that 
suffer stronger collisions, would show more pronounced 
orientational correlations; we study these effects by bin- 
ning the partieles velocities. 

So far we discussed the correlation factor (cos^0)^, 
which is a second moment of the distribution function 
W(cos 0, u, w). Let us now analyze the distribution func- 
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0.032 - • all particles (MD) ' ^ ^ . 

• fast particles (MD) • , 

• slow particles (MD) • . 
0.031- • fast particles (DSMC) 

• J 

0.0 0.2 0.4 0.6 0.8 1.0 

Icos 0| 

FIG. 10; The angular distribution f{\cos6\) for the system 
of rough spheres with et = 0.9, En = 0.9 and q = 2/5 in the 
stationary state. Note that while there is no preferable angle 
between v and <jJ for slow particles, correlations are clearly 
visible for fast particles favouring perpendicular linear and 
angular velocities. 

tion itself. Due to the limited statistics of our numerical 
data we discriminate only between two classes of parti- 
cles: the class of fast particles comprising those particles 
whose linear velocity belongs to the set of the 1/3 largest 
values and whose angular velocity belongs to the set of 
the 1/3 largest values. The class of slow particles is de- 
fined correspondingly as the set of particles whose linear 
velocity belongs to the set of the 1 /3 smallest values and 
the angular velocity belongs to the set of the 1 /3 smallest 
values. In Fig. [TO] we show the distributions /(| cos 6*1) 
for the two classes in comparison with the distribution 
for all particles using both methods, MD and DSMC. In 
both cases we skipped the first 20 collisions per particle 
such that the ratio of temperatures, r, has reached its 
stationary value. For the MD simulation we used a sys- 
tem of TV = 10^ particles at low density (filling factor 
< 1%. Then we averaged over 200 snapshots in distance 
of 1 collision per particle. In case of DSMC we used a sys- 
tem of TV = 10*" particles and made the statistics based 
on a single snapshot. Both results agree very well. The 
angular distribution is almost flat for slow particles and 
cannot be distinguished from the distribution of all parti- 
cles (within statistical accuracy). On the other hand the 
fast particles exhibit a nonuniform distribution with a 
maximum around cos 9 ~ 0. Physically this means that 
the angle 9 between v and uj for slow particles is uni- 
formly distributed within the interval (0, tt), while for 
fast particles it lies preferentially around ~ tt/2. In 
other words, for the particular choice of £t ~ 0.9 and 
En = 0.9 the fast particles tend to behave like sliced ten- 
nis balls, with uj perpendicular to v. 



VI. CONCLUSIONS AND OUTLOOK 

We have analysed in detail the correlations between 
rotational and translational motion in a granular gas of 



frictional particles. Under the assumption of molecular 
chaos and homogeneity we have developed an analyti- 
cal theory which accounts for the correlations (cos^ 0)^ 
in addition to the rotational R{t) and translational T{t) 
temperature. We have also performed large scale DSMC 
simulations as well as event driven simulations to study 
the evolution of a gas of rough spheres and in particular 
the above correlations. 

We observe that the gas of rough particles always 
relaxes to a steady-state with constant correlation 
(cos^6')^ and constant ratio Too — R{t)/T{t). While 
the relaxation of ^cos^ 6') and r to their steady-state val- 
ues happens on the coUisional time scale, the evolution of 
the rotational and translational temperature in the near- 
smooth limit £t — > 1 is a slow process with a vanishingly 
small rate ~77t~(l — £t)<g;l. Physically, this may be 
explained as follows. In the near-smooth limit the cou- 
pling of the rotational modes to the translational ones be- 
comes very weak. The energy of the rotational motion of 
the particles is almost conserved in collisions and the ex- 
change of energy between the translational and rotational 
degrees of freedom becomes very slow. Consequently the 
rotational temperature as well as the translational tem- 
perature have a slowly decaying component, governed by 
this weak exchange of energies. However both temper- 
atures decay with the same slow timescale so that their 
ratio, r, is stationary - after it has reached its steady state 
on the fast time scale of a few collisions. Simultaneously, 
(cos^S)^ relaxes to its steady state with a similar rate 
of the order of a few collisions. We conclude that the 
relaxation of the temperature ratio, r, and the angular 
correlations is rapid, - independent of the strength of 
the coupling (1 — St) as long as it is finite. Furthermore 
the correlations persist up to a vanishingly small rough- 
ness and are absent only for perfectly smooth particles, 

= 1, which makes expansions around the smooth limit 
questionable. 

Our main results concern the correlation between the 
directions of rotational and translational velocity in the 
stationary state: The correlations depend sensitively 
on the values of the coefficients of restitution and the 
moment of inertia; for most of the system parameters 
(cos^ 6*) < 1/3, implying that linear and angular ve- 
locities are preferably orthogonal, like in a sliced ten- 
nis ball. Only for a small part of the parameter space 
(cos^ 9) > 1/3, which means that v and uj are preferably 
parallel like in a rifled bullet; the manifold of vanishing 
correlations (in e„, St space) has seemingly zero measure. 
The correlations are more pronounced for strong devia- 
tions from equipartition. 

Our approach can be extended in several directions. 
In the simulations it is straightforward to use more ad- 
vanced models for the coefficients of restitution as func- 
tions of the impact velocity, e.g. [H, IH, HO, lll|. It 
would also be of interest to study the full one-particle 
distribution. Our results already indicate that more en- 
ergetic particles have stronger correlations, but a sys- 
tematic study has yet to be done. Furthermore, one 
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expects to observe correlations not only in very dilute 
gases, but also in rapidly moving denser systems. Our ap- 
proximate analytical theory is based on the assumption 
of homogeneity and the density only enters into the En- 
skog collision frequency, which sets the time scale. Hence 
our results for the stationary state are independent of 
the density. This cannot hold true in a rapidly moving 
dense system, yet we expect to observe correlations as 
well. These could be analysed in a molecular dynamics 
simulation either for a driven [ZJ or undriven system. 
Finally, the observed correlations may have important 
consequences for the stability theory of dilute granular 
flows: they possibly alter the domain of stability of gran- 
ular system with respect to shear fluctuations — the main 
instability of granular flows of smooth particles. 
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APPENDIX: ANALYTICAL CALCULATIONS 

1. Correlation factor 

We present the details of the analytical calculations, 
leading to the three self-consistent equations (|^^ and 
([M)) for T{t),R{t) and b{t). First, we note that he com- 
putation of b{t) or (cos^ 0)^ is severely hampered by the 
denominator in Eq. Fortunately one can carry out 
the calculations with the auxiliary observable 
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<^)t^^E^''^'^2(cOS0,). 
i—1 



(A.l) 



Its relation to our set of observables can be established 
by essentially the same steps as leading from Eq. to 
Eq. (EH): 



30b(t) 



mm 

qinn?a? 



? 9 

qui a 

mm 



(A.2) 



In the case of vanishing correlations we have (A)^ = 0. 
Positive (negative) values correspond to a preference of 
a parallel (perpendicular) orientation. 

Owing to the assumptions of spatial homogeneity and 
molecular chaos it suffices to consider the phase space of 
only a single pair of particles (without loss of generality 
these shall be labeled 1 and 2). Integrating out the spatial 
degrees of freedom and using the definition of the pair 
correlation function 



N{N- 1) / drs . . .drNQNiri, ■ ■ - Vn) = n .92(7-12), 

A.3) 

with n being the number density of the gas (e.g. [6|) we 



obtain 

{iIi+vA)^^iyN III y (n-i;i2)e(-n-t;i2) 

Vi V2 LOi UJ2 

X [\+h{t)vlLjlP2{cOSei)\ 



(A.4) 



X [l + b[t)vlu}lP2 (cos 6*2) 

X fSi2-l) A, 



where h is an arbitrary but fixed unit vector, v 
— 87rna^g2(2a) and we used the shorthand notations 



f— ) 

\2ttTJ 

f I 



\2nR 



3/2 



3/2 



d Vi exp 



' 2T 

d-'ij.cxp ( 



(A.5) 



In the following we will drop the 6^(t)-term stemming 
from the product of the two one particle distribution 
functions pi since it was assumed to be small and we 
only want to go to first order in b{t). 

The calculation of ^612 — 1^ A is obviously rather in- 
volved and, thus, it needs to be broken up to stay 
tractable. It is convenient to introduce relative integra- 
tion variables 



(A.6) 



V = Vi2/V2 V={vi+V2 

uj = u)i2/V2 ft = {ivi + UJ2) / V2 ■ 



The term {ifi^A)^ can be broken up along two different 
principles. First, one can make the dependence on b{t) 
explicit, that is. 



(zi:+A), = ((&i2-i) a)^°^ 

+ 6(<)((6i2-l) a)^'Vo(62) , (A.7) 
where for any function F 

{Ff^=>^j J J J {n-v,2)ei-n-v,2) F (A.8) 

Vi V2 Uli UJ2 

and 



{Ff^=,y I I I I {n-v,2)ei-n-vi2) 



Vi V2 LUi UJ2 



X [vlLjfP2 (cos6ii) -t^ vlu:iP2 {00362)] F. (A.9) 

In order to be able to exploit some further symmetries it 
is advisable to split up the last average again, 



odd 



(A.IO) 
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where 



3^/2 



{n ■ v)0 (— n • v) 



V V n 



{V-nf + iV-uf 
{v-nf + {v-u:f 



F (A.ll) 



involves only even powers of V, and 

(n • v) 8 (— n • v) 



/r^\odd 3\/2 



V V uj n 

X [{V -n) {v-u}) + {v -u) (v-n) 



-{v-v) (n-u) 



in contrast involves only the odd powers of these quanti- 
ties. 

Independently we can write 



A = Aa- Ab/3 



(A.13) 



where 



Aa = J2 ■ '^')^ and Ab=J2 '"^'^^ ■ (^■^'^) 



F (A. 12) First we address the A^-part. Applying the collision 

rule to yields 



bi2 - l) A^ = • u:f + {S-nf + ^ [{h xS)-v]'' + ^ [(n xS) - Vf -V2{d ■ n) (v ■ O) 



9 ? 
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V2 



V2{5 ■ n) {V ■ u}) - V2{5 ■ u>) [v ■ ui) - V2{8 ■ u}) {V ■ VI) + — {v ■ n) {h X 5) ■ V 

qa 

+ yl(^V ■u){hxS)-v + — {v-u}){nxS)-V+ — {V-n){nxS)-V 

qa qa qa 

- —{6 ■n){nx 6) -v - — {S ■uj){nx 6) -V (A.15) 

qa qa 



and invoking the definition of S, Eq. (j4]), we obtain 
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+ 2 



(0) 



+ 2% ^(n • y)^ (n • ujfy ' + 277j^a^ (^[{h x O) • a;]^^' ' + 2^ ^[(n x v) ■ V]') 



^2\ 



(0) 



2\ 



(0) 



q-^a^ 



l2\ 



(0) 



(2ryt - 1) (r;„ _ r;*) - ^ + 



2 g g 



q 



{{n ■ v) {n ■ Lj) {v ■ u>)) 



(0) 



1 ) ((n-y)(n-a;)(V-a;))(°^ (A.16) 



The terms that vanish by symmetry are already left out at this point. The contributions to ( ( 612 — 1 ) A^ 
exactly the same form. 



have 
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For ^(^12 ^ 1^ one finds the following contributions 



^ - ^ - 47yt j ((t; • a;) (V • - ^ (ry„ - rjt) {{n ■ v) (h ■ V) (n • u;) (n • Jl)) 



odd 



+ (1 - 2r]t) {{fi ■V){h- n) {v ■ io)) 



odd 



([(n X i;) • V] [(n x $7) ■ lj]) 



odd 



((n a;) (V • . (A.17) 



y2 („ X rj)^y°^ 



Correspondingly, the A^-part may be written as 

_ 1) A.y"^ ^ - 1) (2. - If (.^ (. X <V ^ _ 1 

+ ^ (2,. - 1)^ (v^ in X vff + M (y2 x vff + 4,. - 1) (v^u^f^ 
+ 4 {^l - ^„ ~ + ^0 ((n . u'f^ + ^-I^(yJl-l^ {rjl - ^„ - + rj,) ((n • vf [n x f^)^)^"^ 
+ ^ - l) (2.7* - 1) ([(n X .) . n]'f' + ^{vl-Vn- + m) {{n ■ vf (n x vf) 

+ 2if,a^ ((n X 0)2a;2)'"' + 2r^la^ ((n x nf n^^^ + - 1^ ((n x O)^) 



(0) 



(0) 



(0) 



(A.18) 



The contributions to ( [ &12 ^ 1 ) A^ ) again are formally equivalent to the above expression. This leaves us with 



'12 
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477, 



^ (2r7t - 1) ((n • (n • f2) (-y • y))°^^ + ^ ([(n x v) ■ uj] [(n x J^) • V]) 

- ^ - ((n • t;) (n • V) (n • a;) (n • 0))°'^'^ . (A.19) 
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We have now reduced the problem to the tedious but 
straightforward calculation of a considerable number of 
averages. This task is best suited for a computer al- 
gebra system and thus we only tabulate the results. 
To simplify the notation we introduce the abbreviations 
D = v^T/m-K, f = T/m, and R = R/I 
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([(n X t;) • = -2z>f2 (A.20f) 

{{n ■v)ih- uj) (v ■ Lj)f"'' = -2vfR (A.20g) 
((n ■V){n- Lj) {V ■ w))'"^ = -vTR (A.20h) 
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(0) 
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{v^u;^)^°^ = -l2vfR (A.20m) 
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(n X ri)^r22) = -lOzJi?^ 
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(0) 



(n X n)") ' ■ = -8i>^2 



(0) 



(n X n)^ (n X t;)^ ) = -AvTR 
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(A.20r) 
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(A.20u) 

(A.21a) 
(A.21b) 



„ „ \ even _ _ 

(n • vy (n xny) = ei^T^i?^ (A.21j) 



(n X rj)^\ = ivTR^ 



„ even _ ^ 

r2^(nxrj)^) =7vTR^ 



{fi X r2)''V = 8f>™3 



2 2 \ 

(n X d) (n X r2) \ =0 



((^; ■L^)[V ■ rj))"^** = -20Df^R^ (A.22a) 



((n u;) (F • = -lODT^R^ (A.22b) 



((n • V) (n • rj) (t; • = -TDT^R^ (A.22c) 



(A.21k) 
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(A.21m) 

(A.21n) 
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((n V) (n • a;) (n • 0))°'^^ = -Avf^ie 

(A.22d) 



([(n X • V][{n x • u)])"'^'^ = (A.22c) 



((n • w) (n • n) (i; • V'))°'^'^ = ^2DT'R' (A.22f) 



([(n X v) ■ u}][{h X f2) • -buT^R' . (A.22g) 



-QDf^R^ 
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2. The correction terms for the temperatures 

To calculate dT/dt ~ {iL+T)^ one essentially proceeds 
along the same lines of reasoning as detailed above. First 
of all, it is again advantageous to write the corrections to 
the Gaussian distribution function explicitly, that is, 



(0) 



(ii:+r), = ( 612-1 T) 612-1 



(1) 



where 
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f rj^a" (n x $7) 
+ rit {2rit -l)av-{fix Q.) (A.24) 
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and 



The term ^612 — 1^ is already known [2^ and the 
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